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Glacier Ice Mass Fluctuations and Fault Instability in Tectonically Active Southern 
Alaska by Jeanne M. Sauber and Bruce F. Molnia 

Across southern Alaska the northwest directed subduction of the Pacific plate is 
accompanied by accretion of the Yakutat terrane to continental Alaska. This has led to 
high tectonic strain rates and dramatic topographic relief of more than 5000 meters within 
15 km of the Gulf of Alaska coast. The glaciers of this area are extensive and include 
large glaciers undergoing wastage (glacier retreat and thinning) and surges. The large 
glacier ice mass changes perturb the tectonic rate of deformation at a variety of temporal 
and spatial scales. We estimated surface displacements and stresses associated with ice 
mass fluctuations and tectonic loading by examining GPS geodetic observations and 
numerical model predictions. Although the glacial fluctuations perturb the tectonic stress 
field, especially at shallow depths, the largest contribution to ongoing crustal deformation 
is horizontal tectonic strain due to plate convergence. Tectonic forces are thus the 
primary force responsible for major earthquakes. However, for geodetic sites located < 
10-20 km from major ice mass fluctuations, the changes of the solid Earth due to ice 
loading and unloading are an important aspect of interpreting geodetic results. 

The ice changes associated with Bering Glacier’s most recent surge cycle are 
large enough to cause discernible surface displacements. Additionally, ice mass 
fluctuations associated with the surge cycle can modify the short-term seismicity rates in 
a local region. For the thrust faulting environment of the study region a large decrease in 
ice load may cause an increase in seismic rate in a region close to failure whereas ice 
loading may inhibit thrust faulting. 

Prior to the 1979 St. Elias earthquake (M=7.2), the main thrust zone below the 
study region had been locked since the 1899 earthquakes and strain had been 
accumulating. During this same time period ongoing wastage of southern Alaska glaciers 
was 100’s of meters up to almost a 1 km. We used estimates of ice thickness decrease to 
calculate the changes in the fault stability margin around the region of the 1979 St. Elias 
earthquake and aftershocks. Our results suggest that the cumulative decrease in the fault 
stability margin due to ice wastage between 1899 and 1979 was large and would promote 
thrust faulting. Since earthquake hazard evaluations are based on the paleoseismic 
history of the region, for glaciated areas the concurrent glacial history should be 
considered in this evaluation. 
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Abstract 

Across the plate boundary zone in south central Alaska tectonic strain rates are 
high in a region that includes large glaciers undergoing wastage (glacier retreat and 
thinning) and surges. For the coastal region between the Bering and Malaspina Glaciers 
the average ice mass thickness changes between 1995 and 2000 range from 1-5 m/yr. 
These ice changes caused solid Earth displacements in our study region with predicted 
values of -10 mm to 50 mm in the vertical and predicted horizontal displacements of 0 - 
10 mm at variable orientations. Relative to stable North America, observed horizontal 
rates of tectonic deformation range from 10-40 mm/yr to the north-northwest and the 
predicted tectonic uplift rates range from approximately 0 mm/yr near the Gulf of Alaska 
coast to 12 mm/yr further inland. The ice mass changes between 1995-2000 resulted in 
discernible changes in the GPS measured station positions of one site (ISLE) located 
adjacent to the Bagley Ice Valley and at one site, DON, located south of the Bering 
Glacier terminus. In addition to modifying the surface displacements rates, we evaluated 
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the influence ice changes during the Bering glacier surge cycle had on the background 
seismic rate. We found an increase in the number of earthquakes (Ml > 2.5) and 
seismic rate associated with ice thinning and a decrease in the number of earthquakes and 
seismic rate associated with ice thickening. These results support the hypothesis that ice 
mass changes can modulate the background seismic rate. 

During the last century, wastage of the coastal glaciers in the Icy Bay and 
Malaspina region indicates thinning of hundreds of meters and in areas of major retreat, 
maximum losses of ice thickness approaching 1 km. Between the 1899 Yakataga and 
Yakutat earthquakes (M w = 8.1, 8.1) and prior to the 1979 St. Elias earthquake (M s = 
7.2), the plate interface below Icy Bay was locked and tectonic strain accumulated. We 
used estimated ice mass change during the 1899 - 1979 time period to calculate the 
change in the fault stability margin prior to the 1979 St. Elias earthquake. Our results 
suggest that a cumulative decrease in the fault stability margin at seismogenic depths, due 
to ice wastage over 80 years, was large, up to ~2 MPa. Ice wastage would promote thrust 
faulting in events such as the 1979 earthquake and subsequent aftershocks. 
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1. Introduction 


Across southern Alaska the northwest directed subduction of the Pacific plate, v p 
= 53 mm/yr [DeMets et al., 1994], is accompanied by accretion of the Yakutat terrane to 
continental Alaska ( v accr = 33-44 mm/yr) [Plafker et al., 1994]. This has led to high 
tectonic strain rates and dramatic topographic relief of more than 5000 meters within 15 
km of the Gulf of Alaska coast. The glaciers of this area are extensive (Figure 1) and 
many of them have undergone significant advances and retreats of kilometers to tens of 
kilometers on time scales ranging from decades to thousands of years [Molnia, 2003]. 
The repeat times of large earthquakes in this region are hundreds of years [Nishenko and 
Jacob, 1990]. 

In Scandinavia and interior Canada the tectonic deformation rates are low but the 
change in glacier thickness during late Pleistocene deglaciation were large (up to ~3 km) 
over regional to continental spatial scales. During and following late Pleistocene 
deglaciation moderate earthquakes occurred in the stable interior of Canada but large 
earthquakes occurred in Scandinavia (see summary and papers in topical volume Stewart 
et al., 2000). In southern Alaska tectonic deformation rates are high but the recent 
change in ice load is smaller (< 1 km) and local to regional in extent. Our data suggest 
that fluctuation of Alaska glaciers over the last 100 years may influence the timing of 
southern Alaska earthquakes, but tectonic forces largely control the rate and distribution 
of ongoing crustal deformation. 
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In this study we present new estimates of the surface displacements associated 
with ice mass changes following the 1993-1995 Bering glacier surge (1995-2000) and we 
evaluate the influence that associated stress changes may have on the background seismic 
rate. Also, we use an estimate of the cumulative ice thickness decrease between 1899 
and 1979 to calculate predicted stress changes along the main thrust zone that slipped in 
the 1979 St. Elias earthquake (M s = 7.2). 

Although glacial, sedimentary, and tectonic processes in southern Alaska have been 
well studied significant questions remain to be answered [Jaeger et al., 2001]. 
Specifically in this study, we address in Section 2 “How large are recent ice mass 
fluctuations between the Malaspina and Bering Glaciers”? Then in Section 3, “What 
are the tectonic displacement rates and stresses”? Followed by, “How do we quantify the 
influence of ice mass changes on earthquake faulting?” in Section 4. Next, in Section 5 
we address, “What are surface displacement changes and incremental stresses associated 
with recent ice mass fluctuations?” and “Do these ice mass fluctuations alter the seismic 
rate of background seismicity?” Finally, in Section 6 we ask, “What influence could ice 
mass changes have had on large earthquakes in the Icy Bay region?” 

2. Glacier Ice Mass Fluctuations in Alaska 

The last major glaciation episode in Alaska is correlated with late Wisconsin 
fluctuations of the Laurentide Ice Sheet. Radiocarbon ages show onset of glaciation at 24 
Kyr BP, peak glaciation at 18 Kyr BP, deglaciation beginning at about 13.5 Kyr BP, and 
retreat from the continental shelf complete by about 10 Kyr BP [Molnia, 1989; Hamilton, 
1994]. The ICE-4G ice-unloading model includes a simple representation of Alaskan 
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deglaciation [Peltier, 1993, 1994]. The ICE-4G predicted uplift caused by late 
Pleistocene deglaciation was about 2 mm/yr near VYAK (Figure 2). An asthenosphere 
viscosity less than 10 Pa s inferred from post-glacial rebound studies in southern 
Alaska [Sauber et al,. 2000; Larsen et al, 2003; Tamisiea et al., 2003], as well as other 
subduction zones [Ivins and James, 1999; James et al., 2000], is lower than the viscosity 
derived from global postglacial rebound studies. In southern Alaska relaxation times 
would be <100 years. Therefore, the ongoing uplift due to late Pleistocene deglaciation 
is probably smaller than suggested by the ICE-4G model because most of the rebound 
would have occurred much earlier. 

Following late Pleistocene deglaciation, additional large glacier ice fluctuations that 
were either regional or global (Little Ice Age) in extent have been suggested. A major 
regional glacier ice advance, widely documented in Alaska, culminated 2600 to 2800 
years ago throughout the Cordilleran region [Porter and Denton, 1967]. It is interesting 
to note that the age of marine terrace 13 in the Icy Bay region (Figure 1) has a midrange 
calendar age of 2533 yr BP [Plafker et al, 1981; Sauber et al., 2000]. A long period of 
milder climate followed in which glacier volume and extent significantly diminished. The 
next major documented regional advance occurred during the Little Ice Age. 

As summarized by Molnia[2003], during the Little Ice Age the area of glacier 
coverage was at least 10% greater than the ice today. Retreat from the Little Ice Age 
maximum position and thinning began as early as the first half of the eighteenth century 
(pre-1750) in the southeastern Alaska Coast Mountains. In the Chugach-St. Elias 
Mountains region retreat generally began later (post-1850) and today more than 95% of 
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the glaciers extending below 1500 m are currently undergoing thinning, retreat, and or 
stagnation. 

2.1 Surge-Retreat Cycle of Bering Glacier, Alaska in the Last 100 Years 

A significant retreat of the Bering Glacier began sometime between 1905 and 
1910 [Molnia and Post, 1995]. The recession of the Bering Glacier has been interrupted 
by at least six twentieth century surges [Molnia and Post, 1995; Muller and Fleisher, 
1995]. A surging glacier develops a periodically oscillating state characterized by a long 
period of slow build-up of ice in the middle and upper parts of the glacier with retreat in 
the terminus region (intersurge or “quiescent phase” of Bindschadler, 1982), followed by 
a surge [Budd and Mclnnes, 1974; Kamb et al., 1985]. During a surge, ice is transported 
from the upper reaches (the surge “reservoir” region) of the glacier, its surface locally 
lowers by tens of meters as the ice is transported down glacier, where the glacier thickens 
and the terminus frequently advances (the surge “receiving” region). 

During the 1993-1995 Bering Glacier surge, a large transfer of ice from the Bagley 
Ice Valley to the Bering Glacier terminus region occurred. The results of extensive 
efforts to study the 1993-1995 Bering Glacier surge have provided some constraints on 
its timing, spatial extent, and ice thickness changes [ Lingle et al., 1993; Molnia, 1993; 
Krimmel, 1994; Molnia et al., 1994; Molnia and Post, 1995; Sauber et al., 1995, 2000; 
Roush, 1996; Herzfeld and Mayer, 1997; Fatland , 1998; Fatland and Lingle, 1998; 
Molnia, 2003]. Earlier, Sauber et al. [1995, 2000] calculated the elastic displacement of 
the solid Earth due to ice mass changes during the surge and compared them to Global 
Positioning System (GPS) measurements at sites near the Bering Glacier. 
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Since the cessation of terminus advance in mid-September 1995, the Bering Glacier 
returned to actively thinning and retreating. Between 1995 and 2003 the terminus of the 
Bering Glacier lobe has retreated up to 6.5 km and the maximum surface thinning in the 
piedmont lobe exceeded 100 m. Based on repeated aircraft laser profiles along a 
centerline profile, Arendt et al, [2002] estimated an average yearly thinning rate for 
Bering Glacier’s piedmont lobe of 3.1 m/yr between 1995 and 2000. During the latter 
part of the 20 th century, parts of the Bering Glacier's piedmont lobe thinned > 5 m/yr 
[Molnia and Post, 1995 \ Molnia, 2003]. 

Following the surge, in the surge reservoir region, transport of ice from upstream in 
the eastern Bagley Ice Valley has resulted in a thickening of ice near GPS site ISLE 
(Figure 1). Field observations made in August 2000 by B. Molnia suggest that between 
1997 and 2000, in the region south and east of ISLE, the glacier ice thickened by 10-20 
meters while immediately north of the station ISLE, the Jefferies and Tana Glaciers have 
thinned by 1-3 m/yr. The Jefferies and Tana Glaciers were not affected by the Bering 
Glacier surge. 

2.2 Glaciers of Icy Bay and the Malaspina Glacier 

The present Icy Bay is the result of twentieth century retreat of ice mass that filled the 
basin {Molnia, 1977, 1989; Porter, 1989; summary in Molnia, 2003). At the onset of the 
twentieth century an expanded Guyot-Malaspina Glacier that extended several kilometers 
into the Gulf of Alaska filled the entire basin. Retreat began prior to 1910 and continues 
today. During the twentieth century, the glaciers that filled Icy Bay have retreated to 
form an approximately 40 km long, triangular-shaped bay, with four separate fiord arms 
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at its head (Figure 1). Porter [1989] calculated that the mean rate of retreat in Icy Bay 
was ~1 km/yr for the period 1904 to 1926 and ~0.4 km/yr for the period 1926 to 1989. 

Elevations of recent glacial deposits, moraines, and trimlines, show that in many 
recently deglacierized areas of Icy Bay, former ice thicknesses were significantly greater 
than 300 m above m.s.l. (see Table 3, Sauber et al., 2000). Near the present location of 
the terminus of Tyndall Glacier (Figure 1), the 1959 USGS 1:250, 000-scale topographic 
map [USGS, 1983], which was based on 1957 aerial photography, shows > 800 m of 
glacier ice above sea level at a location that presently is ice free. Several hundred meters 
of ice were also lost from below sea level. Since the early twentieth century, when retreat 
began in Icy Bay, the maximum ice thinning in parts of the bay has exceeded a kilometer. 
Although terminus retreat began around 1910, ice wastage in the form of ice thinning 
probably began sooner. 

During most of the twentieth century, the margin of the Malaspina Glacier has been 
stagnating, down wasting, and retreating. Thinning of the piedmont lobe has been 
estimated to be 85 to 129 m in the last 100 years (Sharp, 1961, Lingle et al. 1999; 
Tangbom, 1999;summary in Molnia, 2003). Aircraft laser altimetry measurements by 
Arendt et al. [2002] between 1995 and 2000 indicate an average thinning rate of 
approximately 0.9 m/yr across the Malaspina piedmont lobe. 

3. Tectonic Stress and Strain Orientation and Magnitude 

The ongoing rate of tectonic deformation over the last 20 years has been estimated 
from trilateration, very-long-baseline interferometry (VLBI) and GPS measurements 
[Savage and Lisowski, 1988; Ma et al., 1990; Sauber et al., 1993, 1997, 2000]. The rate 
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and orientation of tectonic stress and strain is due primarily to subduction of the Pacific 
plate and collision of the Yakutat terrane with interior Alaska (Figure 2). Within the 
data uncertainties, the VLBI measurements (1984-1990) at Cape Yakataga (VYAK in 
Figure 2) show a constant velocity for all three components (North, East, and Up) except 
for an offset due to the 1987/1988 Gulf of Alaska earthquakes. In the next section we 
discuss the recent (1993-2001) GPS measurements as well as a summary of our 
numerical modeling of the tectonic component of the surface velocities. 

3.1 Ongoing Tectonic Displacement Rates Across the Study Region 

GPS measurements were made between 1993 and 2000 across a portion of southern 
Alaska (Figure 2). The GPS data were processed using the GAMTT/GLOBK software 
[King and Bock, 1997; Herring, 1997]. Details of observations made in 1993, 1995 and 
1997 are given in Sauber et al. [1997; 2000]. Subsequently, in August 2000, we 
conducted a photo-survey of the Bering Glacier and lower reaches of the Bagley Ice 
Valley to estimate ice thickness change and we made GPS observations at stations ISLE 
and DON (Figure 2). The horizontal velocities for the full time period were estimated 
relative to an ITRF2000 North America fixed reference frame (Table 1, Figure 2). We 
used an approach similar to that discussed by Sauber et al. [2000] to estimate the station 
velocities. To illustrate the evolution in the time series of measurements at DON and 
ISLE, we estimated station position for each observation year as well. 

Relative to stable North America, the measured horizontal rates of tectonic 
deformation across the region given in Figure 2 are 10-40 mm/yr. The predicted tectonic 
uplift rates range from approximately 0 mm/yr near the Gulf of Alaska coast to 12 mm/yr 
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further inland. The observed uplift rates are significantly higher at some stations and in 
part reflect uplift associated with either glacial wastage or the Bering Glacier surge cycle. 

We used a two-dimensional plane strain finite element method (TECTON, Melosh 
and Raefsky, 1981) to model the tectonic deformation rates. The finite element grid 
across this subduction zone plate boundary was modified from Cohen [1996] and 
includes a shallow dipping subducting slab and both an oceanic crust-mantle and a 
continental crust-mantle (see Figure 5 in Sauber et al. 2000). The horizontal velocities 
of sites located > 20 km from major glacier fluctuations are consistent with a model that 
includes the viscoelastic response to a downgoing Pacific plate that is locked at shallow 
depths as well as a component of slip (15 mm/yr) perpendicular to the orientation of plate 
convergence [Sauber et al., 1997, 2002]. Variations in the surface deformation rate due 
to the viscoelastic response to major earthquakes, as well as downdip creep below the 
seismogenic zone, may contribute as well. However, until we obtain more permanent 
GPS sites in this region it is difficult to isolate these contributions. 

3.2 Stress State and Principal Stress Directions 

A number of large earthquakes have occurred in the last 30 years within or near the 
study region given in Figure 1 and 2 (1970, Pamplona zone M w = 6.7; 1979, St. Elias, M s 
= 7.2; 1987-1988, Gulf of Alaska, M s = 6.9, 7.6, 7.6, Figure 2). The static stress drop for 
these events ranges from moderate to high (~2 to 10 MPa) [Hwang and Kanamori, 1992; 
Estabrook et al., 1992; Sauber et al., 1993; Doser et al., 1997]. Although stress levels 
on individual faults are highly variable, these earthquakes suggest that much of the region 
given in Figure 1 is close to failure. 
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The earthquake focal mechanisms, offshore in situ borehole measurements and 
regional geological evidence [Hottman et al., 1979, Estabrook and Jacob, 1991; Plafker 
et al., 1994; Doser et al., 1997] have been used to estimate stress directions. For the 
coastal region between Icy Bay and Cordova (CRDC in Figure 2), a horizontal north- 
south to northwest-southeast orientation is suggested for the maximum stress (aO; east of 
Icy Bay Oi may rotate to more northeast-southwest orientation. A minimum stress, G 3 , 
that is vertical has been assumed for the dominantly thrust faulting environment. Strike- 
slip faulting associated with a vertical intermediate stress has been suggested near the 
Contact fault just north of the Bagley Ice Valley [Savage and Lisowski, 1988]. Based 
on borehole failure observed in offshore wells between Icy Bay and Kayak Island 
Hottman et al. [1979] estimated the approximate magnitude of principal stresses. Their 
findings for depths between 2700 m and 4000 m were extrapolated to obtain estimates for 
the maximum CTi of about 65 MPa and a minimum 03 of about 18 MPa at a depth of 5 km 
[Sauber et al., 2000]. 

4. Shear Failure and Fault Stability Criterion 

The state of stress in southern Alaska is considered to be the result of the 
superposition of the plate tectonic stresses, the stress due to overlying rocks, the pore 
fluid pressure, and the stress due to ice mass changes. In this section, we briefly discuss 
the fault failure criterion used to evaluate the influence of stress changes due to ice mass 
fluctuations. 

4.1 Mohr-Coulomb Failure Criterion 

The Mohr-Coulomb criterion for brittle shear failure in rock is described by 
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t = T 0 + pc„ (1) 


where t is the shear stress necessary to induce failure on a fault plane, t c is the inherent 
shear strength of the fault, p is the coefficient of friction of the fault surface, and a n is 
the normal stress on the fault [Jaeger, 1969]. The locus of shear t 0 and normal a n stress 
components on a suite of faults of various orientations is a Mohr circle whose center is 
the average between the maximum Oi and minimum <33 principal stresses and whose 
radius is (<7i - a 3 )/2 (Figure 3). 

4.2 Fault Stability Margin 

Following the approach given in Wu and Hasegawa [1996], we estimated the 
fault stability margin and the change in that quantity due to ice mass changes (Figure 3). 
For a region in which the Mohr circle representing the stress state in southern Alaska is 
close to, but not in contact with the Coulomb failure envelope, we have the situation 
given in Figure 3. Here we assume that thrust faulting is the primary mode of crustal 
faulting as indicated above in Section 3. Since we are interested in the effect of glacial 
ice mass changes on faulting potential, we calculated the changes from an initial state t Q 
until some time after glacial change, t. The shortest perpendicular distance from the 
failure envelope to the outer Mohr circle is defined at the Fault Stability Margin (FSM, in 
Figure 3) and it is related to the principal compressive stresses by: 

FSM = |3 [p (oi + a 3 ) + 2r 0 ] - V 2 (01 - a 3 ) (2) 

(from Wu and Hasegawa 1996) where |3 = [arc tan (p)]/ 2p. 

The variation of FSM (dFSM) will indicate whether ice mass loading will tend to 
trigger failure or promote stability: 
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dFSM = 1/2 { [a (t 0 ) - a 3 (t o)] ■ [oi (t ) - c 3 (t )]} + 


(3a) 


mid (t ) + o 3 (t )] - [0! (t o) + a 3 (t o)]} (3b) 

Changes in the deviatoric stress or changes in the mean stress will move the 
Mohr circle. The first term (3a) represents the change in the mean stress associated with 
the stress change; the second term indicates the change in the deviatoric stress. Since we 
calculated the difference between an initial state and a later state, static pressure terms 
that affect all the principal stresses in the same way and do not vary in time, are 
automatically subtracted out. Thus, dFSM is independent of pore fluid pressure and the 
overburden stress for the lithostatic case [Wu and Hasegawa, 1996]. In the sections 
below we summarize our estimates of the surface displacements and stresses due to ice 
mass changes and we evaluated their influence on fault failure. 

5. Observed and Predicted Surface Displacement Rates and the Seismic Rate 
Associated with the Bering Glacier Surge Cycle 

The spatial (tens of kilometers) and temporal (1-2 years) scales, as well as the 
magnitude, of loading and unloading associated with the Bering Glacier surge cycle is 
similar to water reservoir impoundment and subsequent water level fluctuations. For 
loads of this spatial scale the response of the Earth is assumed to be primarily elastic or 
poroelastic (see the summaries by Simpson et al., 1988; Scholz, 1990; Talwani, 1997). 
Since the stress increase from an elastic load drops off rapidly with distance, seismicity 
fluctuations in these cases are concentrated in the immediate vicinity of the load changes, 
and earthquake sizes tend to be small, since the stress gradients are high. In the following 
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section we show the geodetic time series for two sites close to the Bering Glacier and we 
present the predicted contribution that ice fluctuations make to the observed station 
positions. Next, we examine the seismicity and seismic rate for portions of the 
intersurge, surge, and post-surge time periods. Here we tested the hypothesis that ice 
unloading promotes thrust-faulting earthquakes, whereas, ice loading inhibits thrust- 
faulting earthquakes. 

5.1 Observed and Predicted Surface Displacements 

As summarized in Section 2 a glacier that periodically surges undergoes a long period 
of slow build-up of ice in the middle and upper parts of the glacier with retreat in the 
terminus region. This is referred to as the intersurge in Figures 4-6. At some point this 
phase is followed by a surge where ice is rapidly transported from the surge reservoir 
region to the surge receiving region [Budd and Mclnnes, 1974; Bindschadler, 1982, 
Kambetal., 1985]. 

The Bering Glacier’s surge reservoir and receiving areas each have a geodetic site 
that was observed during the 1993 to 2000 time period. In Figure 4 the time series of 
station position is given for the North, East, and Up components for DON and ISLE. We 
assumed that the tectonic deformation rate was constant between 1993 and 2000. Earlier 
geodetic measurements at sites located > 20 km from areas of major ice mass fluctuations 
support this assumption [Ma et al., 1990; Sauber et al„ 1993, 2000]. The departure from 
a constant rate of change, especially evident in the vertical component at DON and the 
east component of ISLE, is attributed to ice loading or unloading. We used the ice 
thickness changes between 1995 and 2000 summarized in Section 2 to calculate the 


14 


predicted elastic displacement (Figure 5) following the approach given in Farrell [1972] 
and Sauber et al. [1995]. 

Station DON is located south of the Bering Glacier terminus and was observed in 
1993, 1995, 1997, and 2000 (Figure 4a). The vertical component most dramatically 
demonstrates subsidence associated with ice loading of the ice receiving area in Bering 
Glacier’s piedmont lobe, rapid uplift due to ice thinning and terminus retreat following 
the end of the surge (1995-1997), and then moderate uplift associated with continued ice 
thinning and terminus retreat (1997-2000). Due to an estimated, average thinning rate of 
3.1 m/yr across the ablation zone of the Bering Glacier between 1995 and 2000 [Arendt et 
al., 2002], station DON is predicted to move to the southeast (Figure 5a) and up (Figure 
5b). The predicted tectonic uplift rate is only 0-5 mm/yr so the ice loading is especially 
evident; whereas the predicted horizontal tectonic rate is quite high (almost 40 mm/yr to 
the north-northwest). 

Station ISLE is located near the ice reservoir region in the Bagley Ice Valley and it 
was observed in 1993, 1995, and 2000 (Figure 4b). During the surge in 1993-1995, ice 
thinning in the adjacent Bagley Ice Valley caused ISLE to be displaced toward the 
northeast and up [Sauber et al., 1995, 2000], Following the surge, some time between 
1995 and 1998, the ice began to thicken in the Bagley Ice Valley but most of the Bering 
Glacier continued to thin. The combination of ice loading south of ISLE and the slower 
glacial wastage ice of surrounding glaciers caused ISLE to be displaced down slightly 
and to the east-southeast relative to the 1995 position (Figure 5). The predicted tectonic, 
vertical component is larger (about 12 mm/yr) than DON, whereas the predicted 
horizontal component at ISLE is smaller than DON, about 20 mm/yr to the north- 
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northwest. As shown in Figure 2, the 1995-2000 horizontal velocity of ISLE is indeed 
more north directed than the surrounding station velocities which are oriented to the 
north-northwest. 

5.2 Seismicity and Seismic Rate 

In addition to surface displacement associated with ice mass changes, we hypothesize 
that with ice unloading thrust earthquakes would be promoted and with ice loading 
earthquakes would be inhibited. We plotted earthquakes of magnitude Ml > 2.5 between 
1973-2003 (Figure 6). Included on the seismicity plots are three reference areas: (1) the 
surge reservoir region in Bagiev Ice Valley near GPS site ISLE, (2) the lower portion of 
the Bering Glacier and the adjacent coastal plain that includes GPS site DON, and (3) the 
Cape Yakataga region where in the last 100 years only minor scale glacier fluctuations 
have occurred within 20 km of the VYAK. The Cape Yakataga region showed very few 
earthquakes throughout the 1973-2003 time period with the exception of a couple of 
earthquakes around the time of the 1979 St. Elias earthquake and one during the surge 
time interval (Figure 6). 

For the Bagley Ice Valley (surge reservoir region), we examined the seismic rate for 4 
time periods. For Jan. 1973 - Dec. 1992, following the 1967 surge and before the 1993 
surge, thus “intersurge”; for Jan. 1993 - Dec. 1995, the “surge” time period; for Jan. 1996 
- Dec. 1998, the” post-surge” readjustment period, and finally ice build-up in the 
reservoir region, “intersurge”. In addition to examining the change in the number of 
earthquakes associated with stress changes [Stein, 1999], we calculated the cumulative 
seismic moment and the seismic moment rate because this quantity would be less 
sensitive to the detection threshold. Most of the earthquakes (Ml > 2.5) in the reservoir 
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region between 1973 and 2003 occurred during the surge time period or shortly 
afterwards in the post-surge time period (Figure 6). In addition to an increase in the 
number of earthquakes during the surge and post-surge time period, relative to the longer 
intersurge (1973-1992) time period, there was an order of magnitude increase in the 
seismic moment rate during the surge (Table 2). 

Cumulative ice thinning in the reservoir region, during the surge, was on average ~50 
meters across a 50 km region. For a short time following the surge, ice draw-down 
continued. The maximum dFSM associated with ice thinning during the surge was about 
-0.5 MPa; the maximum change in the stress field would be at the surface in a localized 
region just below the ice thinning. However, at some time between 1996 and 1998 the 
ice began to thicken. From April 1998 to August 2003 no earthquakes (Ml > 2.5) 
occurred in the surge reservoir box. 

For the Bering Glacier terminus region (surge receiving area), there were few events. 
One event occurred during the early intersurge time period and one Ml = 4.3 occurred 
during the later intersurge time period. During the intersurge time periods the Bering 
Glacier piedmont lobe undergoes terminus retreat as well as ice thinning. The one event 
during the surge time period occurred very early in April 1993 before the Bering Glacier 
terminus advanced during the surge. In the surge receiving area ice thickening and 
terminus advance (1993-1995) added an average of 50-100 meters of ice; the dFSM 
would be up to 1.0 MPa and thrust earthquakes would be inhibited. Then, following five 
years of ice thinning and retreat after the surge, there was one event in 2001 (Ml = 4.3). 

6. Ice Mass Fluctuations and the Occurrence of Large Earthquakes 
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Prior great earthquakes in the region given in Figure 1 occurred on September 4, 1899 
(M w = 8.1, Yakataga) and September 10, 1899 (M w = 8.1, Yakutat) [. Plafker and Thatcher, 
1982; Nishenko and Jacob, 1990]. These events were large enough that they may have 
relieved much of the strain accumulation in the region referred to as the Yakataga seismic 
gap [Nishenko and Jacob, 1990]. The first event apparently was an interplate earthquake, 
which involved faulting along the Chugach-St. Elias fault zone. It produced coastal uplift 
of 1-3 m between Icy Cape (western mouth of Icy Bay) and Cape Yakataga (VYAK, in 
Figure 2) [Plafker and Thatcher, 1982; Jacoby and Ulan, 1983], The second event 
caused massive uplift along the northwest shore of Yakutat Bay [Tarr and Martin, 1912]. 

The St. Elias earthquake of February 28, 1979 (M s = 7.2), occurred on a down-dip 
extension of the Chugach-St. Elias fault system (Figures 1 and 2). Detailed modeling of 
seismic body and surface waves suggest a mainshock location to the northwest of Icy Bay 
(Figure 1) and an epicentral depth of 24 km. A best-fitting source mechanism indicates 
initial underthrusting on a shallow north-dipping plane that later changed to a more 
steeply northeast dipping plane with a right-lateral component [Estabrook et al., 1992]. 
The ongoing tectonic strain accumulation across the region is shown in Figure 2 and 
discussed in Section 3. Next, we discuss the predicted viscoelastic response of the Earth 
to ice mass changes preceding the 1979 earthquake. 

6.1 Changes in the Fault Stability Margin Prior to the 1979 St. Elias Earthquake 

Since retreat of the Icy Bay glaciers started around the beginning of the twentieth 
century, we calculated the change in the fault stability margin associated with estimates 
of ice mass changes in the -80 years prior to the 1979 earthquake. Earlier, Sauber et al. 
[2000] used a scaled, idealized profile of ice thickness change to calculate predicted 
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surface displacements rates for comparison to a nearby geodetic site, AMBR (Figure 2). 
The results from that study suggest a preferred viscosity of 5x 10 19 Pa s for depths > 38 
km. Here we estimated ice thickness changes based on the data summarized in Section 2. 
We estimated ice thickness changes along a profile extending from the western mouth of 
Icy Bay to the center of the aftershock zone (Figure 6). The estimated ice changes along 
this profile were 300-800 meters in the coastal region and 100 meters southeast of the 
center of the aftershock zone. The average ice change value over a 10 km region was 
input as surface (node) force to the same viscoelastic finite element model used to 
estimate the tectonic contribution to surface deformation rates. The predicted surface 
displacements north of the maximum ice thickness change were upward and away from 
the trench (northwest, north, and northeast); whereas, just south of the region of 
maximum ice change, surface displacements were directed upward and toward the 
southwest, south, and southeast (Figure 5). 

The change in the fault stability margin (dFSM) was calculated on the surface and at 
depth along the main thrust zone. We calculated dFSM in two different ways. In the 
first case we assumed an initial stress state of zero following the 1899 earthquakes. In 
the second case we assumed the stress magnitudes given for a depth of 5 km (discussed in 
Section 3) as an initial state. As expected the dFSM was the same in both cases. The 
dFSM due to the elastic response to ice mass changes was the largest component in the 
viscoelastic response. The dFSM due to total ice changes over ~80 years were up -3.6 
MPa at the surface and -2.2 MPa along the main thrust zone. The largest change occurs 
in the region just north of the maximum ice thickness change where the maximum 
increase to Gi as well as a decrease in <53, the vertical stress, occurs. These changes 
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increase the size of Mohr circle, bringing the region closer to failure (Figure 3). The 
location of dense earthquake occurrence associated with the 1979 event (Figure 6) is 
north of the region of large ice thickness change. 

7. Discussion and Summary 

The large glacier ice mass changes associated with glacial wastage and surges 
perturb the tectonic rate of deformation at a variety of temporal and spatial scales. We 
estimated surface displacements and stresses associated with ice mass fluctuations and 
tectonic loading by examining GPS geodetic observations and numerical model 
predictions. Although the glacial fluctuations perturb the tectonic stress field, especially 
at shallow depths, the largest contribution to ongoing crustal deformation is horizontal 
tectonic strain due to plate convergence. Tectonic forces are thus the primary force 
responsible for major earthquakes. However, for geodetic sites located < 10-20 km from 
major ice mass fluctuations, the changes of the solid Earth due to ice loading and 
unloading are an. important aspect of interpreting geodetic results. 

The ice changes associated with Bering Glacier’s most recent surge cycle are 
large enough to cause discernible surface displacements. Additionally, ice mass 
fluctuations associated with the surge cycle can modify the short-term seismicity rates in 
a local region. For the thrust faulting environment of the study region a large decrease in 
ice load may cause an increase in seismic rate in a region close to failure whereas ice 
loading may inhibit thrust faulting. We did not observe a change in the seismicity rate 
associated with smaller seasonal ice mass fluctuations (i.e. late summer versus late 
winter). 
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Prior to the 1979 St. Elias earthquake, the main thrust zone below Icy Bay had 
been locked since the 1899 earthquakes and strain had been accumulating. During this 
same time period ongoing wastage of southern Alaska glaciers was 100’s of meters up to 
almost a 1 km. We used estimates of ice thickness decrease to calculate the changes in 
the fault stability margin around the region of the 1979 St. Elias earthquake and 
aftershocks. Our results suggest that the cumulative decrease in the fault stability margin 
due to ice wastage between 1899 and 1979 was large and would promote thrust faulting. 
Since earthquake hazard evaluations are based on the paleoseismic history of the region, 
for glaciated areas the concurrent glacial history should be considered in this evaluation. 
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Table 1 . Site Velocities in a North American Fixed ITRF 2000 Reference Frame 
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Table 2. Earthquakes and Seismic Moment Rate as a Function of Time in the Bering 
Glacier’s Surge Reservoir Region 


Date 

#Eqs 

Seismic Moment 

(dyne -cm) 

Seismic Moment Rate 
(dyne-cm/yr) 

Jan. 1973-Dec. 1992 

2 

2.4 x 10 21 

1.2 x 10 20 

Jan. 1993-Dec. 1995 

7 

2.2 x 10 22 

7.3 x 10 21 

Jan. 1996-Dec. 1998 

3 

3.5 x 10 20 

1.1 x 10 20 

Jan. 1999- Aug. 2003 

0 

0 

0 


Figure Captions 


Figure 1. The August 9, 2003 true color, Moderate Resolution Imaging 
Spectroradiometer (MODIS) image from NASA’s Aqua satellite of the southern Alaska 
coastal region. The MODIS image shows the Malaspina piedmont lobe in the east and 
the Bering piedmont lobe in the west. The solid star shows the location of the 1979 St. 
Elias epicenter. The image is courtesy of the MODIS Rapid Response Project at 
NASA/GSFC. The inset in the lower left hand comer shows the region of Alaska 
covered by Figure 1. 

Figure 2. Location of geodetic sites occupied between 1993 and 2000 in which the 
horizontal velocities have been estimated (solid circle with site abbreviations from Table 
1). The vector velocities of sites relative to an ITRF2000 North America fixed reference 
frame are given with error ellipses representing regions of 95% confidence. The shaded 
regions (with year) indicate the rupture zones of major earthquakes in southern Alaska 
during the twentieth century. The NUVEL-1A predicted, horizontal velocity of the 
Pacific plate relative to a fixed North American plate (52 mm/yr at N14°W, DeMets et 
al., 1994) is shown by the vector in the lower right hand comer. Faults (solid lines): BRF 
= Border Ranges, CF = Contact, CSEF = Chugach-St. Elias, DF = Denali, FF = 
Fairweather, PFZ = Pamplona fault zone, TF = Totschunda, TFZ = Transition, *MW = 
Mount Wrangell, the youngest volcano in the western Wrangell volcanic zone. The 
dotted line indicates the region given in Figure 1. 



Figure 3. Mohr circle, Coulomb failure line, Fault Stability Margin (FSM) and thrust 
faulting mode (modified from Wu and Hasegawa, 1996). The stress state is defined by 
the principal (compressive) stress components: the maximum principal stress Ci, the 
intermediate principal stress (a 2 ), and a minimum effective stress 03 . For the thrust 
faulting mode considered here, the maximum (oO and intermediate principal stress (a 2 ) 
are horizontal and the least principal stress is vertical (a 3 ). 

Figure 4. North, east, and vertical station position evolution for the station DON (a) and 
ISLE (b). The station positions for both sites are predicted to move to the northwest and 
up due to tectonic forcing. For the 1993-1995 surge interval, the predicted influence of 
ice unloading near ISLE and ice loading near DON are given by Figure 2 and Plate 1 in 
Sauber et al. [2000]. For the 1995-2000 post-surge time interval, the predicted influence 
of ice changes from section 2 of this paper is given in Figure 5a, b. 

Figure 5. Predicted displacement of the solid Earth due to ice changes between 1995 
and 2000 represented by disk loads, (a) Predicted horizontal elastic displacement field 
(millimeters) of the solid Earth due to ice thinning and retreat in the ablation zone of the 
Bering Glacier (red disks) and ice loading in the Bagley Ice Valley (blue disks). The 
value in meters from the bottom left to the top right: - 20 , - 20 , - 20 , - 20 , - 20 , - 20 , - 20 , - 20 , 
-20, -20, -20, -20, -20, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -15, -10, -10, - 
10, -10, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, 10, 20, 20, 20, 20 , 20 , 
20, 20, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5. (b) A contour plot of the predicted uplift and 


subsidence (millimeters) of the solid Earth associated with the loading/unloading given in 
Figure 5 a. 

Figure 6. Earthquakes of M L > 2.5 between 59° N and 61° N and 139° W and 144° W 
between 1973 and 2003. The long-dashed rectangles indicate the surge reservoir in the 
Bagley Ice Valley and the surge receiving area of the Bering Glacier terminus; the short 
dashed line indicates the Cape Yakataga rectangle. Fault abbreviations are as given in 
Figure 2. GPS stations are given by open squares. 
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